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INTRODUCTION 


The  purpose  of  this  document  is  to  describe  the  statistical  distribution 
models  that  have  been  adopted  for  use  in  a  model  evaluation  demonstration 
software  package  developed  for  the  Acoustic  Model  Evaluation  Committee  (AMEC). 
Most  of  the  statistical  tests  and  measures  of  accuracy  that  have  been 
incorporated  into  this  package  are  documented  elsewhere  (see  references  in 
McGirr  [1979]).  However,  as  an  aid  in  assessing  the  statistical  nature  of 
empirical  data,  a  means  has  been  provided  to  allow  graphical  comparison 
against  curves  generated  by  theoretical  distribution  models.  For  example,  the 
empirical  probability  density  function  (pdf)  of  data  derived  from  propagation 
loss  measurements  can  be  compared  graphically  against  curves  generated  by  the 
log-gamma  pdf.  As  another  example,  residual  errors  obtained  by  taking 
differences  between  propagation  loss  measurements  and  predictions  can  be 
compared  graphically  against  curves  generated  by  a  scaled  Student's  t  pdf. 
Also,  Kolmogorov-Smirnov  tests  have  been  included  in  the  software  package  to 
allow  testing  empirical  cumulative  distribution  functions  (CDF)  against  either 
the  normal  (or  Gaussian)  CDF  or  the  log-gamma  CDF.  Certainly  the  popularity 
enjoyed  by  both  the  normal  and  the  Student's  t  distributions  precludes  the 
necessity  of  elaborate  discussion  concerning  either  one.  The  log-gamma 
distribution,  on  the  other  hand,  does  not  accede  to  general  application  as  do 
the  normal  and  the  Student's  t.  What  is  more,  the  log  transformation  adopted 
here  is  reversed  in  sign  from  that  typically  cited  in  the  acoustics 
literature,  and  consequently  expanded  discussion  is  desirable. 

Section  2  discusses  the  impetus  for  choosing  the  gamma  distribution  vice 
either  the  exponential  distribution  or  the  "sum-of-exponentials"  distribution. 
Maximum  likelihood  parameter  estimates  are  discussed  as  an  alternative  to 
estimates  based  on  the  method  of  moments.  Some  basic  properties  of  the  gamma 


distribution  are  reviewed,  and  one  of  the  parameters  of  the  gamma  distribution 
is  given  an  interpretation  based  on  a  comparison  with  the  distribution  of  the 
sum  of  two  exponentially  distributed  variables. 

Section  3  discusses  the  log-gamma  distribution  that  is  obtained  by 
applying  a  log-transformation  to  gamma  distributed  intensities.  The 
particular  transformation  applied  here  is  reversed  in  sign  from  that  usually 
reported  in  the  acoustics  literature,  and  an  explanation  for  this  reversal  is 
provided.  Sasic  properties  of  the  log-gamma  density  and  distribution 
functions  are  presented,  and  a  brief  discussion  of  parameter  estimation  is 
also  included. 

The  report  concludes  with  a  few  remarks  pertaining  to  useful  extensions 
of  the  AMEC  software  package,  and  possible  applications  of  the  log-gamma  and 
other  distribution  functions  to  the  problem  of  placing  bounds  on  propagation 
loss  predictions. 
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2.  THE  GAMMA  MODEL 


Statistical  models  of  intensity  fluctuations  for  signals  (and  noise) 
propagating  through  an  ocean  medium  have  been  presented  and  analyzed  by  Dyer 
[1970,  1973].  More  recent  articles  by  Mikhalevsky  and  Dyer  [1978]  and  by 
Mikhalevsky  [1979]  have  expanded  the  distribution  theory  pertaining  to  signal 
and  noise  statistics  under  various  conditions.  Of  particular  interest  here 
are  probability  density  and  distribution  functions  that  are  useful  in 
characterizing  the  statistical  nature  of  propagation  data. 

For  harmonic  signals  under  well-developed  multipathing,  Dyer  [1970] 
determines  the  probability  density  function  (pdf)  of  intensity  (short-time 
average  of  mean  square  pressure)  to  be  exponential  (Rayleigh  amplitude). 
Letting  x  denote  intensity,  then  x  has  pdf 

p(x)  =  A  1  exp(-x/A),  x>0,  \>0,  2.1 

2  2 

with  mean  px  =  \  and  variance  a x  =  \  .  Coming  from  a  radically  different 

point  of  view,  Dozier  and  Tappert  [1978]  develop  a  statistical  theory  of 

normal  mode  amplitudes  in  a  random  ocean.  For  acoustic  frequencies  less  than 

about  250  Hz  and  under  saturation  conditions,  they  arrive  at  the  same  pdf  in 

the  limit  as  the  number  of  modes  increases  without  bound. 

When  the  multipath  intensities  are  not  identically  distributed  but 

can  be  sorted  into  M  groups,  then  each  x  ,  m  =  1,2,...,M,  is  exponentially 

m  M 

distributed  with  mean  The  pdf  of  the  sum  x  =  I  xm,  say  f(x),  takes  a  form 
resembling  a  weighted  sum  of  exponentials  [Dyer,  1973].  That  is, 

f(x)  =  *  (wm/yexp(-xAm)  2.2 

m-1 

M  M  2 

with  mean  g  =  2  \  ,  variance  o  c  =  l  A  ,  and  "weighting  factors"  given  by 
*  m=l  m  *  m=l  m 
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w'^ 


M 

n  (i- 

i/m 


AA>- 


On  the  other  hand,  if  the  Am  =  A  for  all  m,  then  the  sum  has  the  gamma  pdf, 
say 


g(x) 


x 


M-l  -x/A 
e 


AMr(M) 


? 

with  mean  u  =  MA  and  variance  o 
Mx  x 

Graybi 1 1  [1963]). 


2.3 


MA 


(see  for  example  p  126  of  Mood  and 


Allowing  for  distinct  group  means  has  intuitive  appeal,  and  when  the 
underlying  assumptions  prevail  the  pdf  f(x)  is  probably  a  more  realistic  model 
than  the  gamma  pdf.  There  is  an  unappealing  aspect  of  this  model,  however, 
which  stems  from  the  less  than  satisfactory  means  available  to  estimate  the 
set  of  Am  from  data.  Ambiguous  estimates  of  the  Am  may  be  obtained  from 


Thus,  unless  there  exists  a  priori  information  identifying  M-l  of  the  Am, 
additional  information  must  be  acquired  from  higher-order  sample  moments. 
Procedures  based  on  higher-order  moments  are  computationally  prohibitive  for 
M>2  due  to  the  range  of  values  typical  of  intensities. 

For  the  special  case  of  only  two  groups  (M  =  2)  the  method  of  moments 
becomes  less  impractical,  but  such  a  method,  though  simple  conceptually,  is 
not  without  drawbacks.  Consistent  estimates  of  A^  and  A^  may  be  obtained  from 
(see  p  270  of  Lindgren  [1978]) 

Aj  +  A^  =  x 

2  2  2 
and  Aj  +  A2  =  S 
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T IT 


_  2 

where  x  and  S  are  the  sample  mean  and  variance. 
If  A^  >  A 2>  say,  then 


A1  =  *s(x  +  ^2S2  -  x2) 
and  A^  =  ^(x  -  ^2S2  -  x2) 

Since  the  A^  are  real  and  positive,  the  sample  variance  is  necessarily 
constrained  to  lie  in  the  interval  (x2/2,  x2).  The  likelihood  that  this 
constraint  could  be  easily  violated  suggests  that  the  "sum  of  exponentials" 
model  lacks  the  robustness  necessary  to  tolerate  possible  discrepancies 
between  the  model's  underlying  assumptions  and  reality. 

2. 1  The  Gamma  as  a  Generalized  Exponential 

Even  if  a  reasonable  technique  to  estimate  the  Am  could  be  devised,  there 
is  still  no  clear  cut  way  to  determine  M,  except  in  "controlled"  experiments. 
Thus  there  is  ample  incentive  to  consider  the  gamma  pdf  expressed  in  the  form 


g(x) 


»v*VxA 

\vr(») 


2.4 


as  an  appropriate  density.  Note  that  M  has  been  replaced  by  v.  The  purpose 
of  this  exchange  is  to  emphasize  that  this  parameter  is  no  longer  restricted 
to  integer  values.  An  interpretation  of  v  is  given  in  Section  2.4. 

The  idea  of  using  the  gamma  pdf  as  a  generalized  exponential  is  not  new. 
In  his  development  of  the  statistics  of  energy  fluctuations,  Rice  [1954] 
suggests  the  gamma  pdf  expressed  in  the  form 

p(E)  =  an+1Ene'aE/r(n+l) 
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as  an  appropriate  density  of 


One  of  the  advantages  of  the  gamma  distribution  is  that  estimates  of  \ 
and  v  may  be  obtained  from  first-order  statistics.  For  a  sample  of  size  N  the 
likelihood  function  takes  the  form 


L(A,v)  =  n  x  v_1e'Nx/A/(AVr(v))N, 
n=l  n 


where  x  is  the  arithmetic  sample  mean 


x  =  VT  1  x  . 
Nn=l  n 


If  x  denotes  geometric  sample  mean,  that  is 


X  -  (nx/'», 


£nL  =  N[(v-l)£nx  -  xA  -  v£n\  -  £nf(v)]. 


*  / 


Maximum  likelihood  estimates  of  X  and  v  are  obtained  by  solving  the  likelihood 
equations 


3£nL 

3X 

3£nl_ 

3v 


=  N(x/X2  -  v/X)  =  0 
=  N(£nx  -  £nX  -  4i(v))  =  0, 


where  i|>(v)  =  d  £nr(v)/dv  is  the  digamma  (or  psi)  function.  Thus  the  pair  of 
equations 


and 


vX  =  x 


£nv  -  ip( v )  =  £n(x/x) 


2.5 


2.6 


provide  an  iterative  scheme  to  obtain  estimates  of  both  X  and  v  from 
first-order  sample  statistics. 

This  method  of  parameter  estimation  offers  a  computational  advantage  over 
methods  using  second-order  moments.  The  only  moment  that  need  be  calculated 
in  intensity  space  is  x.  The  geometric  mean  is  actually  calculated  in  dB 
space  as  an  arithmetic  mean.  That  is,  under  the  transformation  y  = 
a  £n(x/vX),  where  a  =  -10  log  e, 

£nx  =  £nvA  +  |a|  ^  i  I  y 
N  ■'n 

=  £nx  +  I  a  | -1  J  I  yp 

and  hence 

£nv  -  iJj(v)  =  I  a |  1  jqlyn  =  y/|a|.  2.7 

At  first  glance  this  expression  does  not  appear  to  offer  a  straight¬ 
forward  approach  to  the  problem  of  estimating  v.  However,  since  the  left  hand 
side  of  this  equation  is  independent  of  sample  data,  it  may  be  plotted  against 
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v,  thus  providing  a  simple  graphical  procedure.  This  procedure  is 
demonstrated  in  section  3.5. 

2.2.2  Method  of  Moments 

An  alternative  procedure  ifj.  offered  by  the  so-called  "method  of  moments" 

i 

(see  p  278  of  Lindgren  [1978]),  j wherein  px  and  a^2  are  expressed  in  terms  of 
the  parameters  v  and  X  and  then  {equated  to  the  corresponding  sample 
statistics.  Thus, 

Px  =  vX  =  x, 
and 

ox2  =  vX2  =  S2 , 

from  which 

v  =  x2/S2 .  2.8 

Rice  [1954]  suggests  similar  estimates,  except  instead  of  using  the  usual 
sample  statistics,  x  and  S2,  he  expresses  the  moments  in  terms  of  the  spectral 
density  function  of  the  process.  More  recently  Frisk  [1978]  demonstrates  the 
use  of  this  method  in  estimating  a  parameter  analogous  to  v  for  the  so-called 
M-distribution  (which  may  be  obtained  from  the  gamma  by  making  the 
transformation  z  =  x/v). 

The  maximum  likelihood  esti  mation  procedure  is  generally  preferred  to 
other  methods  available.  Lindgren  notes,  however,  that  the  method  of  moments 
usually  yields  "efficient"  estimates  that  serve  as  useful  first-order 
approximations. 

2.3  Some  Properties  of  the  Gammji  Distribution 
2.3.1  Gamma  Density  Curve 

If  X  =  1  the  gamma  pdf  assuHnes  the  simple  form 


8 


1 


v-1  -X 

g(x)  =  - — — ,  x>0,  v>0. 


A  further  simplification  obtains  when  v  is  confined  to  the  set  of  positive 


integers,  say  v  =  n,  n  =  1,2,...,  so  that 


n-1  -x 

g  (X)  =  * - ^  X>Q)  n=i ,2,3,...  . 

°  (n-1)! 


Curves  of  this  simplified  form  are  illustrated  in  fig  1  for  n  =  1,  2,  3,  and 
4.  For  n  =  1  the  curve  starts  with  a  slope  of  -1,  for  n  =  2,  it  starts  with  a 
slope  of  +1,  a  for  n  >  2  all  curves  have  initial  slope  zero.  The 
exponential  curve  (n  =  1)  has  no  interior  points  of  inflection,  the  n  =  2 
curve  has  an  inflection  point  at  x  =  2,  and  for  n  >  2  each  curve  has  two 
points  of  inflection  at  x  =  n  -  1  ±V n  -  1.  The  modes  (or  most  probable 
values)  occur  at  x  =  n  -  1.  When  \  is  not  unity,  the  expected  value  of  x  is 
nA  and  the  mode  occurs  one  A-unit  to  the  left  at  (n  -  1)A. 


n  =  1 


x"'1  e'x 


Figure  1.  Gamma  pdf  curves  for  n  =  1,2,3,  and  4. 


*  ✓ 


2.3.2  Moments 


The  moments  of  gamma  variates  may  be  obtained  from 


oo  n+v-1  -x/A. 

E  x"  =  ;  2 - s - 32 

o  \  r(v) 

Making  the  change  of  variable  z  =  x/A  yields 


2. 10 


Exn  =  [An/r(v)]  X  zn+v_1e"zdz  =  Anr(n+v)/r(v) 


2.11 


(=MX) 


Thus  the  first  two  moments  are 
Ex  =  vA 
and 

Ex2  =  v(v+l)A2 

from  which  the  variance  o^2  is 


2.12 


2.13 


ox2  =  E(x  -  gx)2  =  Ex2  -  E2x  =  vA2. 


2.14 


Higher-order  central  moments  may  be  obtained  from 


E(x  -  u  )"  =  ; 

o 


oo  .  ,n  v-1  -x/A  . 

n  _  /-  ( x  —  u  )  x  e  dx. 


Avr(v) 


Making  the  same  change  of  variable  as  above  and  simplifying  yields 


2.15 


■  mx>"  =  rTv)  ki0  (nk)Ak(-Mx)n'kr(v  +  k). 


2.16 


The  third  central  moment,  say  m^,  is  then 


m3  =  2vA  , 

so  that  the  coefficient  of  skewness,  a3  is  given  by 
a3  =  m3/m23/2  =  2v‘1/2. 


2.17 


2.18 
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2.19 


r 


Similarly  the  fourth  central  moment  is 

4 

m^  =  3v(v  +  2)A 

and  hence  the  coefficient  of  kurtosis,  a^,  is  given  by 
a4  =  m4/m22  =  3  ^“7^ 


2.20 


Note  that  as  v  gets  large  these  coefficients  tend  toward  Gaussian  values,  that 
i  s 

1  im  (Yj  =  0  2.21 

V  -»  oo 

and 

lim  a4  =  3  2.22 

v  “>  03 

The  rates  of  approach  to  these  limits  are  demonstrated  in  fig  2. 


2.3.3  Gamma  CDF 

The  gamma  CDF  may  be  expressed  as 


Gv(x)  =  f 
0 


t 


y-ytA 

Avr(v) 


dt 


2.23 
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which  for  arbitrary  v  is  not  generally  expressable  in  closed  form.  Tables  of 
the  incomplete  gamma  function  [Pearson,  1934]  may  be  used,  however,  making  use 
of  relationships  and  tables  in  Abramowitz  and  Stegan  (A&S)*[1965]  may  be  more 
convenient.  For  integer  v,  say  v  =  n,  a  finite-term  series  obtains,  that  is 

Gn(x)  =  1  -  e^x/^e'*7*,  2.24 

where 

n  k 

e  (x)  =  I  xVk!  (A&S ,  p  262). 
n  k=0 

When  n/2  is  an  integer,  tables  of  the  chi-square  distribution  may  be  used. 
Making  the  transformation  z  =  2x/\  yields  the  pdf 


p(z) 


z 


n-1  -z/2 
e 

2nT(n) 


so  that  replacing  n  by  k/2  and  z  by  x2  yields 


2.25 


p(x2) 


xk-2e~X2/2 

2k/2r(k/2) 


2.26 


the  chi-square  pdf  with  k  "degrees  of  freedom."  Thus, 

Gn(x)  =  px2(2x/\|2n)  =  1  -  QxZ(2x/\|2n)  2.27 

where  P  2  and  Q  2  are  defined  on  p  940  of  A&S,  and  tabulations  of  Q  2  begin 
XX  X 

on  p  978. 

2.3.4  Remarks  Concerning  Numerical  Evaluation 

When  v  does  not  take  an  integer  or  half-integer  value,  then  either  tables 
must  be  interpolated  or  algorithms  must  be  prepared  for  computer 


*  Subsequent  citings  to  this  reference  are  abbreviated  A&S. 
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implementation.  The  algorithms  used  in  the  AMEC  demonstration  package  are 
based  on  series  and  asymptotic  expansions  that  appear  on  p  262-3  of  A&S.  Now 


x  .v-1  -t/A  x/A  .v-1  -t 

G  (x)  =  J  -  at  =  J  - - —  dt  =  P(v,  x/A) 

o  A  T(v)  o  T(v) 


2.28 


where  P(v,z)  =  y(v,z)/T(v)  =  1  -  r(v,x)/r(v).  The  functions  P(v,z),  y(v,z) 
and  P(v,x)  represent  different  ways  of  expressing  the  incomplete  gamma 
function,  and  are  defined  on  p  260  of  A&S.  Expanding  e  and  integrating 
yields  the  convergent  series  (A&S,  262) 


P(v,z)  =  — -  [1  *  —  —  +  —  -  ...  ]• 

vP(v)  v+1  1!  v+2  2! 


2.29 


Integrating  by  parts  (u  =  e  t,  dv  =  tv  ^dt)  yields  the  slightly  faster  series 


v  -  z  2 

P(  v ,  z )  =  (1  +  —  +  - - -  +  ...],  2.30 

vr(v)  v+l  (v+l)(v+2) 

but  one  that  is  still  slow  for  large  z.  The  asymptotic  expansion 

» 

P(v,z)  =  1  -  *V  1<? -  [1  +  —  +  +  ...  ]  (A&S,  p.  263)  2.31 

r(v)  z  z2 

obtained  by  integrating  by  parts  (u  =  tv  dv  =  e  ^dt)  is  divergent  for  all  z 
and  noninteger  v,  but  for  small  v,  say  v  <  3,  and  "large"  z,  say  z  >  4A,  only 
a  few  terms  are  required  for  sufficient  accuracy. 

2.4  Relationship  to  Sum  of  Two  Exponentials* 

The  purpose  of  this  section  is  to  provide  an  interpretation  of  v  as  it 
ranges  through  values  between  two  successive  integers.  For  propagation  data  v 


*  This  idea  was  suggested  to  me  by  Mr  LK  Arndt. 
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is  expected  to  cluster  about  v  =  1,  or  v  =  2,  corresponding  to  one  or  two 
primary  groups  of  multi  paths.  When  there  are  two  groups  and  the  energy  is 
proportioned  equally  among  them,  their  mean  values  are  equal  and  v  =  2 
exactly.  However,  when  their  mean  values  differ  then  v  takes  on  some  value 
between  1  and  2.  To  see  if  v  can  be  related  in  some  way  to  the  proportioning 
of  energy  among  two  groups  of  multi  paths,  it  is  compared  with  a  parameter  r 
introduced  into  the  "sum  of  exponentials"  pdf. 

Let  x^  and  x2  represent  the  random  variables  associated  with  two  groups 
of  multipaths,  then  each  has  a  pdf  of  the  form 

P^x)  =  A."1  exp(-x/A.),  x  >0,  i  =  1,2,  2.32 

where  each  x^  has  mean  A,..  Without  loss  of  generality  let  A^  >  and 
introduce  the  parameter  r  (0<r<l)  such  that  with  A^  =  A  then  A?  =  rA.  The  pdf 
of  x  =  x1  +  x2  is  given  by  (for  example  see  p  317  of  Parzen  [I960]) 


x 

P(x)  =  /  p1(t)p2(x-t)dt 
0 

=  exp(-x/A)  -  exp(-x/rA)  33 

A(l-r) 

and  E(x)  =  A3  +  A2  =  (l+r)A. 

For  gamma  distributed  variates,  E(x)  =  vA.  Evidently  then  the  parameter 
v  corresponds  to  1  +  r,  at  least  to  first  order.  This  correspondence  does  not 
hold  up  when  higher-order  moments  are  compared,  and  exact  correspondence 
obtains  only  at  the  limit  points  of  r.  Thus  the  implied  linear  relationship  v 
=  r  +  1  is  only  approximate.  Nevertheless,  the  interpretation  is  clear.  For 
the  case  being  considered,  as  v  increases  from  1  to  2  the  partitioning  of 
energy  varies  from  complete  imbalance  (only  one  group)  to  complete  balance 
(two  groups  with  equal  means). 


The  comparison  at  hand  is  illustrated  schematically  in  figure  3. 


Multipath  Partitioning 


One  _ Unequal _ _  Two 

Group  Means  ""Groups 


0  — -  r  - -  1 

e_x/A-e_x/rX 


Figure  3.  Correspondence  between  v  and  r. 


To  see  how  well  the  gamma  pdf  approximates  the  sum  of  two  exponentials, 
let  \  =  1  and  r  =  1/2  so  that  the  sum-of-exponentials  pdf  is 

p(x)  =  e  X~e  2X  =  2e"x(l  -  e"X)  2.34 

1/2 

Under  the  assumed  correspondence,  the  gamma  pdf  has  \  =  1  and  v  =  3/2  so  that 
g(x)  =  x1/2e"x/T(3/2)  =  (2/^t)  x1/2e'x  2.35 

These  curves  are  illustrated  in  figure  4.  The  gamma  pdf  reaches  a  maximum  (x 
=  .5)  slightly  ahead  of  the  other  curve  (x  =  .69),  but  the  curves  tend  to 
coincide  for  x  >  2.5.  Indeed,  at  x  =  10  E(x)  =  15  the  curves  differ  by  only 
about  7  x  10'7.  Thus,  although  this  discussion  is  necessarily  brief  and  deals 
only  with  specifics,  it  demonstrates  the  basis  for  employing  the  simple  gamma 
model  with  v  allowed  to  assume  any  positive  real  value  as  determined  from  a 
sample. 
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X 


2 


3 


iparisan  of  gamma  pdf  with 
i-of-2  exponentials. 
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THE  LOG-GAMMA  DISTRIBUTION 


3. 1  The  Log  Transformation 

The  range  of  intensity  values  typically  encountered  in  ocean  acoustics 
tends  to  be  discouragingly  prohibitive  both  computationally  and  graphically. 
Log  transforming  to  decibels  (dB)  circumvents  possible  computational  problems 
(eg,  underflow)  and  certainly  results  in  more  convenient  graphics  (eg,  60 
units  vice  106).  Thus,  for  example,  propagation  loss  versus  range  plots 
typically  display  loss  in  dB  as  illustrated  in  figure  5. 


Figure  5.  Plot  format  typical  of  propagation  loss  models. 


The  expression  of  acoustic  quantities  in  dB  is  a  well-established 
practice,  so  well  that  there  is  a  natural  inclination  to  transform 
mathematical  "operations"  from  linear  space  to  dB  space  as  well.  As  Dyer 
[1970]  points  out,  however,  the  average  of  log- transformed  data  does  not  equal 
the  log  of  data  averaged  in  linear  space.  For  example,  averaging  in  dB  space 
vice  averaging  in  linear  space  and  then  transforming  can  produce  a  2.5  dB 
error.  The  reason  for  this  discrepancy  is  made  clear  in  section  3.3. 
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Proceeding  to  the  transformation,  let  x  denote  intensity  (actually  the 
ratio  of  intensity  at  a  given  point  to  some  reference  intensity)  then  the  loss 
H  expressed  in  decibels  is  obtained  by 
H  =  -10  log  x. 

For  a  sample  of  intensities  xn,  n  =  1,2,...,N,  obtained  at  ranges  rn  the 
corresponding  sample  loss  values  are 


Hn  =  H(rn)  -  -101ogxn,  n  =  1,2 _ N. 


If  the  range  interval  spanned  by  the  rn  is  not  too  large  then  "fluctuations 
about  the  mean"  may  be  examined  by  calculating  moments  of  Hn  -  <H>,  where  <H> 

=  -lOlogx  and  x  is  the  arithmetic  mean  of  the  intensities. 

Assuming  that  the  xn  are  independent  random  variables  all  deriving  from  a 
given  population  distribution,  the  corresponding  distribution  in  dB  space  may 
be  obtained  through  the  transformation  H  =  -10  log  x.  The  purpose  of  the 
minus  sign  is  clarified  in  the  next  section.  Mathematical  manipulations  are 
somewhat  less  cumbersome  when  the  natural  log  is  substituted  for  the  common 
log.  Hence  the  transformation  is  expressed  as  H  =  a  £nx,  where  a  =  -10  log  e. 
Since  the  expected  value  of  gamma  distributed  variates  is  vX,  the 
transformation  of  interest  here  takes  the  form 

y  =  a  £ n(x/vX).  3. 1 

Thus  the  gamma  pdf 

v-l-x/A 

9(x)  =  ~  r  - -  ,  x  >  0,  3.2 

a  r(v) 

transforms  to 

f(y)  =  g(x(y))/|dy/dx| 


|a|r(v) 


{exp  [y/a  -  exp(y/a)]}v  , 

-oo<y<oo. 


3.3 
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Despite  its  foreboding  form,  this  pdf  is  not  as  exotic  as  its  appearance  might 
suggest.  In  fact  it  is  routinely  applied  in  the  analysis  of  extreme-value 
problems. 

Interestingly,  if  the  motive  for  making  a  transformation  were  based  on 
the  desire  to  obtain  a  variate  with  constant  variance,  this  same  pdf  (or  one 
of  similar  form)  would  result.  That  is,  suppose  a  transformation  z  =  h(x)  is 
sought  such  that  z  has  (approximately)  constant  variance.  Expand  h(x)  about  x 
=  p  and  assume  that  to  a  good  approximation 

h(x)  =  h(p)  +  h‘ (p)(x  -  p). 

Then  E[h(x)  -  h(p)]2  =  [h'(p)]2  E(x  -  p)2  , 
or  oz2  =  [h'(p)]2ox2, 

and  if  ox2  can  be  expressed  as  some  function  of  p,  say  V(p),  then  setting 
~l/2  ~ 

h'(p)  =  [V(p)]  yields  oz2  =  1.  Hence  the  desired  transformation  is 
provided  by 

x 

z  =  h(x)  =  J  dp/VvTil). 

The  variance  of  gamma  distributed  variates  is  v\2,  so  that  for  constant  v 

h(x)  «  /  ~  a.  £nx, 

which  implies  that  the  transformation  yielding  a  variable  with  constant 
variance  is  logarithmic.  A  discussion  of  this  technique  (referred  to  as  the 
"angular  transformation")  may  be  found  in  Brownlee  [1965]. 

3.2  Intensity- Loss  Orientation 

The  material  in  this  section  is  intended  to  orient  the  reader  to  the 
special  application  of  the  log-gamma  distribution  to  the  analysis  of 
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propagation  loss  data  as  generally  viewed  by  those  engaged  in  propagation  loss 
modeling.  The  majority  of  propagation  loss  models  under  consideration  by  the 
AMEC  for  evaluation  offer  an  output  option  that  generates  a  plot  (or  plots)  of 
loss  versus  range.  These  plots  are  usually  oriented  like  the  one  presented  in 
figure  5,  that  is,  with  increasing  values  of  loss  directed  downward.  Since 
the  gamma  pdf  is  highly  skewed,  its  orientation  relative  to  such  plots  is 
important. 

Figure  6  displays  the  orientation  of  distributions  in  linear  (or 
intensity)  space  and  in  dB  space,  under  the  transformation  y  =  a£n(x/vX)  where 
a  =  -10  log  e.  Values  of  x  falling  below  vX  (the  mean)  yield  positive  values 
of  y  (high  loss)  corresponding  to  the  "right"  tail  of  f(y).  Similarly,  values 
of  x  above  vX  yield  negative  values  of  y  (low  loss)  corresponding  to  the 
"left"  tail  of  f(y).  Thus  the  random  variable  y  is  oriented  in  the  same  sense 
as  propagation  loss. 

In  terms  of  sample  data, 

yn  =  Hn-<H>,  3.4 

where  <H>  =  -10  log  x  and  x  denotes  the  arithmetic  mean  of  intensity  data  (x  = 
2x  /N).  Let  y  denote  the  mean  calculated  in  dB  space,  then 

*  =  fan  =  SIHn  -  <H>  =  S  -  <H>-  3.5 
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3.3  Moments 


The  gamma  pdf 


v-1  -x/A 

g(x)  =  -  ,  x  >  0 

a  r(v) 

under  the  transformation  y  =  a  £nj'x/vA)  takes  the  form 
v  ,  y/a  i 

-  (eJ  )t  ,  -oo<y«». 

lalT(v) 

The  first  moment  of  y  is  given  by 


f(y)  = 


Ey  =  X  yf(y)dy 


or,  with  the  prospect  of  more  tractable  mathematics  in  x  space, 


Ey  =  Ex[a  £n(x/vA)] 

00  I 

=  [a/AvT(v)]  J  £n(x/vA)xvt1e‘xA  dx 


Making  the  change  of  varfabie  z  ~jx/A  gives 
“  -l  I  - 

Ey  =  [a/F(v)]  f  £n(z/v)  zv  At(>  zdz 
o 

=  [a/F(v)]  /  (£nz  -  £nv)  ,J:V  Ae  zdz 


=  [a/r(v)]  {  /  £nz  zv‘V;|dz  -  £nv  X  zv"Vzdz}. 
0  0 

oo  , 

Now  r(v)  =  /  zv  Ae  zdz, 
o 

and*  =  J  £nz  zv  zdz, 

dv  o 


*  I  thank  Mr  LK  Arndt  for  bringing  this  fact  to  my  attention. 


3.6 


3.7 


3.8 

3.9 

3.10 
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and  hence 


Ey  =  a[r’(v)/r(v)  -  £nv].  3.11 

Since  r'(v)/T(v)  =  d£nl~(v)/dv  =  tj>(v),  the  digamma  (or  psi)  function  (see  eq 
6.3.1,  A&S),  then  the  expression  for  the  first  moment  of  y  simplifies  to 

Ey  =  a[4»(v)  -  £nv].  3.12 

Evaluating  this  expression  at  v  =  1  yields  (see  eq  6.3.2,  A&S) 

Ey  =  a  4»(1)  =  (-4.  343)(-.  5772) 

=2.5  (dB)  3.13 

which  agrees  with  the  result  obtained  by  Dyer  [1970]  for  the  "log-exponential" 
pdf. 

Let  p  =  E(y)  so  that  the  higher-order  central  moments  are  given  by 


00 

E(y  -  p)n  =  r_1(v)  ;  [a  £n(z/v)  -  p]nzv-1e‘zdz,  3.14 

o 

after  making  the  same  change  of  variable  as  above.  Applying  the  binomial 
expansion  to  the  quantity  in  brackets  and  integrating  over  the  sum  yields 


E(y  -  M)n 


V  (J)  <-q)k  jr  £n"  kz  zV  le  *dz 

k=0  o  f(v) 


♦  (-q)n 


3.15 


where  q  =  £nv  +  p/a. 

Expanding  this  expression  for  n  =  2,  3,  and  4  gives 

E(y  -  P)2  =  a2[r"(v)/r(v)  -  2qr'(v)/r(v)  +  q2],  3.16 

E(y  -  p)3  =  a3[r'  "(v)/r(v)  -  3qr' '(v)/r(v)  +  3q2r'(v)/r(v)  -  q3],  3.17 

and 

E(y  -  p)4  =  a4[r(4)(v)/r(v)  -  4qr"  1  (v)/F(v)  +  6q2r"  (v)/r(v) 

-4q3r'(v)/r(v)  +  q4].  3.18 


23 


The  derivatives  nn'(v)  can  be  expressed  in  terms  of  polygamma  functions.  A 
fairly  thorough  account  of  these  functions  is  given  in  A&S  (p  258-260)  with 
tables  for  4^n^(v),  n  =  0,  1,  2,  and  3  beginning  on  p  267.  The  polygamma 
functions  are  defined  as 

*(n)(v)  =  dn  4»(v)/dvn  3.19 


where 


tj>Cv)  =  d  2n  T(v)/dv  3.20 

so  that 

r'(v)  =  T(v)»|)(v).  3.21 

In  the  following  expressions  for  the  derivatives  P^n^(v)  the  v  dependence  is 
suppressed.  From  f  =  nj>, 

r"(v)  =  dr'/dv  =  r>  +  r*'  =  r(4)2  +  <j/')  3.22 

r'"(v)  =  dr"/dv  =  r'(4>2  +  4*')  +  r(2iw>'  +  4>“)  3.23 

=  r(4<3  +  3«w>‘  +  4>‘ ') 

r(4)(v)  =  dr1  "/dv  =  r'(<i<3  +  34^'  +  4/' ')  +  r(3t|>fy  +  3(4.'  )2  +  3441 '  +  ^ ' ') 

=  r(4)4  +  64j24»'  +  44rP*  '  +  3(4>')2  +  41'")  3.24 

Substituting  into  the  central  moment  expressions  given  above  and 
collecting  terms  yields 

E(y  -  m)2  =  a2  V(v)  3.25 

E(y  -  p)3  =  a3  «l>"(v)  3.26 


E(y  -  p)4  =  a4  [3(4»'  (v))2  -*■  4»'"(v)]. 


The  second  central  moment,  or  variance,  for  v  =  1  is  simply 


*  / 


3.28 


E(y  -  M)2  5  °y2  =  a2  *,(l). 

From  eq  6.4.2  and  23.2.24  of  A&S 

+'(1)  =  £(2)  «  ,  3.29 

so  that 

oy  =  (4. 343)(1. 283)  =  5.57,  3.30 

in  agreement  with  well-known  results  (see  Dyer  [1970]). 

The  coefficients  of  skewness  and  kurtosis,  a3  and  or^,  are  given  by 

a3  =  -*,,(v)/[*1<v)]3/2  3.31 

and 

“4  =  3  +  ij<"  1  (v)/[ii)' (v)]2  3.32 

Asymptotic  expressions  for  «|»' ,  4»' 1  and  i|i'  1  ’  are  given  on  p  260  of  A&S  and 
are  (including  the  first  two  terms  only) 

4,1  (v)  ~  v  +  2^  +  •••  (Eq  6‘412’  MS) 

^"(v)  ~  -  -  Js  -  ...  (Eq  6.4.13,  A&S) 

<l»'"(v)  ~  ^  ^  +  ...  (Eq  6.4.14,  A&S) 

Thus  to  lowest  order  in  v 
-1/? 

a3(v)  ~v  -  0  as  v  ■>  ®  3.33 

and 

a4(v)  ~  3  +  v  ^  -»  3  as  v  -»  ®,  3.34 

and  hence  these  coefficients  are  asymptotically  Gaussian. 

3.4  Some  Properties  of  the  Log-Gamma  Distribution 
3.4.1  Density  Curves 

The  log-gamma  pdf  illustrated  in  figure  7  for  v  =  .05(.25)2.5,  exhibits  a 

pronounced  skew,  as  does  the  gamma  pdf.  Under  the  negative  log  transformation 
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V.  ■* 


*  * 


applied  here,  the  skew  is  to  the  right  of  the  mean.  As  v  increases  the 
skewness  becomes  much  less  pronounced.  A  feature  resulting  from  transforming 
about  Ex  =  vA  is  the  invariance  of  the  mode  with  respect  to  v.  That  is,  let 
h(y)  =  v(y/a  -  e^a)  and  let  Cy  =  vv/|a|r(v),  then 

f(y)  =  cveh(y) 
and  hence 

f '  (y)  =  (v/a)(l  -  ey//a)Cy  eh^ 

Setting  f'(y)  to  zero  yields  the  most  probable  value  of  y,  the  mode,  as 

y  =  a  £n  1  =  0 
for  all  v. 

The  second  derivative  takes  the  form 

f ' ' (y)  =  a'2[-vey/a  +  v2(l  -  ey/a)2]Cveh(y) 
and  when  set  equal  to  zero  yields  a  pair  of  inflection  points  for  each  v,  say 

y+  =  a  £n[l  +  1  *  ^4v  *  X].  3.35 

2v 
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As  v  increases  the  distance  between  the  corresponding  pair  of  inflection 
points  decreases.  Thus  the  distance  between  inflection  points  could  serve  as 
a  measure  of  peakedness. 

3.4.2  Log-Gamma  CDF 

Given  that  x  has  the  distribution  Gy(x)  (eq  2.28),  then  under  the 
transformation  y  =  -|a|£n(x/vA)  the  distribution  of  y,  say  F  (y  ),  is  given  by 

Fv(yo^  =  Pr(y  -  y<P  =  Pr("|aUn(x/vA)  <  yQ) 

‘yo/|a| 

=  Pr(x  >  vAe  ) 

-yo/|al 

=  1  -  Pr(x  <  vAe  ) 

-y0/|a| 

=  1  -  Gv(vAe  0  ) 

“yo/,a| 

=  1  -  P(v,  ve  0  ),  3.36 

where  P(v,x)  is  discussed  in  section  2.3.3. 

Plots  of  Fy(y)  are  presented  in  figure  8  for  v  =  . 75( . 25)2. 25. 
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3. 5  Parameter  Estimation 

By  taking  the  transformation  about  the  expected  mean  intensity,  E(x)  = 
vX,  only  the  parameter  v  needs  to  be  estimated  from  samples  in  dB-space.  The 
parameter  X  is  implicitly  determined  through  the  calculation  of  x  in  intensity 
space.  In  practice  X  need  not  be  known  explicitly.  Since  y  =  a  £n(x/vX)  = 
a  £nx  -  a  £nvX  then  denoting  propagation  loss  measured  at  range  rn  by  H^,  the 
corresponding  sample  point  is  yn  =  Hr  -  <H>  where  <H>  =  -10  log  x  and  x  = 

2xn/N.  The  yn  from  a  specified  range  interval  represent  deviations  of  loss 
from  the  interval  mean. 

Assuming  that  the  yn  are  independent  and  identically  distributed  ( i id)*, 
the  maximum  likelihood  estimator  of  v  can  be  determined  from  (see  section 
2.1.1) 

£nv  -  4»(v)  =  -  y/|a|,  3.37 

where 

y  =  Iyn/N  =  H  -  <H>.  3.38 

Figure  9  presents  y  versus  v,  providing  a  graphical  means  to  estimate  v  given 
y.  The  dashed  curve  gives  a  versus  v  and  can  be  used  as  a  consistency  check 
for  v  obtained  from  y,  since  a2  =  N  X(yn  “  y)  is  an  estimate  of  E(y-p)  =  a 
«!»’  (v). 

The  following  examples  illustrate  the  maximum  likelihood  parameter 
estimation  method.  The  first  example  is  taken  from  a  bottom- limited  data  set 
analyzed  by  Pederson  [1980].  The  sample  size  is  a  rather  awesome  5,398  points, 
which  perhaps  accounts  for  the  quality  of  fit  that  is  evident  in  figure  10. 


*  Obtaining  i id  sample  points  is  evidently  an  ancient  artform  handed  down  from 
one  generation  of  data  reduction  specialists  to  another,  and  involves  closely 
guarded  secret  black  arts  such  as  detrending,  decimation,  and  "proper" 
selection  of  interval  length. 
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The  impressive  sample  size  was  obtained  by  combining  range-detrended  data  sets 
derived  from  several  independent  measurement  events.  The  curve  over-plotted 
on  the  histogram  is  a  log-gamma  density  function  corresponding  to  v  =  1.12. 

The  appropriate  value  of  v  can  be  obtained  graphically  from  figure  9  or  by 
iteratively  solving  equation  3.37  for  v  given  that  y  =  2.22. 

An  example  based  on  model -generated  data  is  demonstrated  by  the  histogram 
of  figure  11.  A  subsample  consisting  of  250  points  generated  by  RAYMODE-X  was 
range-detrended  by  subtracting  the  incoherent  outputs  from  the  coherent 
outputs.  The  sample  mean  is  1.7;  thus  from  figure  9  the  estimated  value  of  v 
is  1.3.  The  closeness  of  fit  is  not  particularly  impressive  near  the  mean, 
but  the  degree  of  coincidence  is  certainly  acceptable  along  the  tails. 
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4.  RECOMMENDATIONS  FOR  FUTURE  WORK 


4. 1  Extensions  to  AMEC  Software 

The  log-gamma  and  Gaussian  "graphical  aids"  are  included  in  the  AMEC 
software  package  largely  as  a  result  of  efforts  in  support  of  a  limited 
propagation  data  assessment  project  recently  conducted  by  Pedersen  [1980]. 
Initial  data  snooping  clearly  indicated  that  most  of  the  collected  data 
derives  from  distributions  exhibiting  skewed  tails.  As  a  consequence  there 
was  no  pressing  need  to  consider  other  distribution  possibilities.  Log-gamma 
and  Gaussian  distributions  do  not  always  provide  adequate  fits,  however,  as 
demonstrated  by  Frisk  [1978]  in  his  comparison  of  exponential,  gamma  and 
log-normal  distributions  against  low-frequency,  long-range  propagation  data. 

At  the  lowest  of  three  frequencies  (9.8,  110,  262  Hz)  all  of  the  distributions 
produced  an  acceptable  fit.  At  110  Hz  only  the  log-normal  with  a  =  5.6  dB  was 
adequate,  and  at  262  Hz  none  was  acceptable. 

Frisk  suggests  that  the  Rice  distribution  may  yield  a  reasonable  fit  to 
data  when  log-normal  and  gamma  distributions  fail.  Derivations  of  the  Rice 
pdf  are  given  in  Rice  [1954],  Davenport  and  Root  [1958],  Whalen  [1971]  and 
01'shevskii  [1967].  The  Rice  pdf  may  be  put  in  the  form 


p(v;a)  =  ,e-1/2<»2  *  *2>I0(a,) 

where  v  and  a  are  ratios  of  the  envelope  and  the  sine  wave  amplitude  to  RMS 
noise.  Urick  [1977]  examines  this  distribution  in  terms  of  a  randomicity 
parameter  T  defined  as 
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When  a  -»  0,  T  -*•  1  and  p(v;a)i  reduces  to  the  Rayleigh  pdf,  and  when  both  v 

and  a  are  large  p(v;a)  exhibits  Gaussian  behavior.  Since  T  ->  0  as  a  -»•  ®  then 

') 

T  ranges  between  0  and  1.  As  an  estimate  of  T,  Urick  uses  the  value  that 

i 

yields  the  closest  agreement  between  the  Rice  CDF  and  the  sample  CDF.  This 
approach  appears  to  be  more  efficjient  (computationally)  than  the  maximum 

i 

likelihood  procedure.  That  is,  rieverting  back  to  the  parameter  "a"  (for 

'.j 

notational  convenience),  the  log-|  1 i kel ihood  function  for  a  sample  of  size  N  is 


’  1  ?  '  ? 

£n  L(a)  =  £n  n  vn  -  ±  I  v^j-  Na  /2  +  I  £n  I  (av  ). 
n=l 

Setting  d  £n  L(a)/da  to  zero  yields 

a  =  ^  1  v  I,  (a  v  )/I  (av  ) 

N  n  1  v  n'  ov  ni 

which  must  be  solved  iteratively-  -a  discouraging  calculation  to  contemplate. 


A  transformation  of  the  form  u  =  jv2  results  in  a  non-central  chi-square 

i 

distribution  (p.  113,  Whalen  [19711])  for  which  extensive  tables  and 
computational  algorithms  have  be^'n  developed  (Harter  and  Owen  [1973]).  For 


AMEC  purposes,  however,  confining  T  (or  equivalently  a  non-centrality 
parameter)  to  a  few  selected  valtfies  appears  to  offer  the  most  efficient 
procedure.  8 

The  Rice  distribution  provides  a  connecting  link  between  the  Rayleigh 
distribution  ("slow"  fluctuation  )  and  the  log-normal  distribution  ("fast" 
fluctuations).  Delineating  the  ..pace  of  fluctuations  vertically  with  a 
"phase-fluctuation  strength"  par  imeter  4>  and  horizontally  with  a  "diffraction" 
parameter  A,  Flatte  et  al  [1977]  suggest  the  Rice  distribution  to  characterize 

l 

amplitude  fluctuations  along  the  fuzzy  border  separating  saturated  (*2A>>1) 
and  unsaturated  (4>2A<<1)  regions/  Thus  as  a  "transition  curve"  the  Rice 
distribution  affords  an  alternat  ve  to  "arbitrary"  Pearson  and  Edgeworth 


f 


frequency  curves  (see,  for  example,  Elderton  and  Johnson  [1969]).  Part  of  the 
impetus  for  characterizing  amplitude  fluctuations  with  "standard" 
distributions  vice  "best-fit"  frequency  curves  stems  from  the  desire  to  have  a 
"physics-based"  framework  for  calculating  bounds  on  model  predictions. 

4.2  Prediction  Bounds 

Serious  attempts  to  place  confidence  bounds  on  propagation  loss 
predictions  have  been  avoided  because  the  dependence  of  fluctuations  on  model 
input  parameters  is  not  well  enough  understood.  Moreover,  the  amount  of 
information  contained  in  a  given  set  of  model  inputs  might  not  be  sufficient. 
Eclipsing  such  pessimism,  however,  the  assumption  is  made  that  RMS 
fluctuations  can  be  related  in  some  way  to  model  input  parameters. 

Three  approaches  to  this  problem  are  considered  here.  The  first 
approach,  and  the  most  expedient,  assumes  that  the  deviation  of  propagation 
loss  from  its  mean  value  follows  a  N(0,o2)  law  with  o  =  5.6  dB.  Using  a 
second-pass  scheme,  the  predicted  loss  versus  range  curve  is  divided  into 
range  intervals  according  to  predominant  path  type.  For  each  such  interval 
the  90%  confidence  limits  are  given  by  Hn  ±  1.96a/VN,  where  N  is  obtained  by 
dividing  the  decorrelation  length  into  the  range  interval  length.  A  slight 
variation  to  this  approach  entails  substituting  a  more  "appropriate" 
distribution  for  the  normal  law,  such  as  the  log-exponential,  the  log-gamma, 
or  the  log-Rice  distribution. 

The  other  two  approaches  are  by  comparison  much  more  extensive  in  scope 
in  that  each  requires  considerable  data  analysis  effort.  In  one  of  these 
approaches,  the  dependence  of  RMS  fluctuations  on  environmental-acoustic  (EVA) 
parameters  is  determined  using  regression  methods.  That  is,  the  RMS 
fluctuation  or  at  range  r  is  regressed  on  selected  EVA  parameters  (or 
functions  of  EVA  parameters). 
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The  third  (and  most  elaborate)  approach  applies  Baysian  methods  (see  for 
example  p  234-7  of  Lindgren  [1978]).  Let  y  denote  the  deviation  of  predicted 
loss  H  from  its  mean  <H>  over  a  specified  range  interval,  and  let  w  denote  the 
deviation  of  measured  loss  H  from  its  mean  <H>.  Suppose  that  w  is  found  to 
have  a  distribution  with  pdf  p(w|0)  conditioned  on  a  known  value  of  the 
parameter  0.  The  parameter  0  can  take  the  form  of  an  n-tuple,  although  in 
practice  it  typically  is  either  a  single  parameter  or  a  pair  of  parameters. 

For  example  if  p(w|0)  is  log-exponential  then  0  =  A,  or  if  p(w|0)  is  log-gamma 
then  0  =  (v,A),  or  if  p(w|0)  is  Gaussian  then  0  =  (p,o2).  During  a  given 
measurement  exercise  information  about  0  is  gained,  yielding  a  prior  pdf  for 
0,  say  g(0).  A  posterior  density  function  h(0|w)  is  then  given  by 

h(0|w)  =  g(0)p(w| 0)/p(w) 

where  p(w)  =  J  g(0)p(w|0)d0.  Finally  the  predictive  density  function  f(y|w) 
is  determined  from 

f(y|w)  =  f  f(y |0)h(0| w)d0 

=  f  f(y|0)p(w|8)q(0)d0 
f  p(w|0)g(0)d0 

The  pdf  f(ylw)  along  with  its  estimated  parameters  can  then  be  used  in 
calculating  confidence  limits.  Note  that  both  the  regression  approach  and  the 
Baysian  approach  tend  to  be  "EVA-specific"  since  in  each  case  the  final 
results  are  tied  to  a  particular  set  of  environmental-acoustic  conditions. 
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